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The overlap Dirac operator at nonzero quark chemical potential involves the computation of the 
sign function of a non-Hermitian matrix. In this talk we present an iterative method, first proposed 
by us in Ref. [1], which allows for an efficient computation of the operator, even on large lattices. 
The starting point is a Krylov subspace approximation, based on the Arnoldi algorithm, for the 
evaluation of a generic matrix function. The efficiency of this method is spoiled when the matrix 
has eigenvalues close to a function discontinuity. To cure this, a small number of critical eigen- 
vectors are added to the Krylov subspace, and two different deflation schemes are proposed in this 
augmented subspace. The ensuing method is then applied to the sign function of the overlap Dirac 
operator, for two different lattice sizes. The sign function has a discontinuity along the imaginary 
axis, and the numerical results show how deflation dramatically improves the efficiency of the 
method. 
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1. The overlap operator and the sign function at nonzero quark chemical potential 

The overlap Dirac operator [2, 3] provides an exact solution of the Ginsparg-Wilson relation 
and hence implements chiral symmetry in lattice QCD even at finite lattice spacing. At zero quark 
chemical potential the overlap operator requires the computation of the sign function of the Hermi- 
tian Wilson-Dirac operator, for which efficient methods have been developed [4, 5] . 

To describe QCD at nonzero baryon density (see Ref. [6] for a review), a quark chemical 
potential /i is introduced in the QCD Lagrangian. The massless overlap Dirac operator at nonzero 
jJL was defined in Ref. [7] as 

Dov(Ai) = l + 75Sgn(//w(At)) (1-1) 
with Hvi{lJL) = YsDv^du.). D^{ji.) is the Wilson-Dirac operator at nonzero chemical potential [8] 

[DMU = Sn,n - £ (1 + 7y)f/«,;5„+;,,„ - K- £ (1 - 7,-)^;, /n-y> (1-2) 

7=1 7=1 

where K= 1/(8 + 2m^ ) with negative Wilson mass m ^ G ( — 2 , 0) , 7v with v = 1 , . . . , 4 are the Dirac 
gamma matrices in Euclidean space, 75 = 71727374, and U„,v are SU{3) matrices. The exponential 
factors implement the quark chemical potential on the lattice. For jJ. y^O the argument H^{ijl) 
of the sign function in Eq.(l.l) becomes non-Hermitian, and one is faced with the problem of 
defining and computing the sign function of a non-Hermitian matrix. 

Consider a given matrix A of dimension and a generic function /. Let F be a collection of 
closed contours in C such that / is analytic inside and on F and such that F encloses the spectrum 
of A. Then the function f{A) of the matrix A can be defined by [9] 

f{A) = ^J^f{z){zI-A)-'dz, (1.3) 

where the integral is defined component-wise and / denotes the identity matrix. From this definition 
it is easy to derive a spectral function definition. If the matrix A is diagonalizable, i.e., A = UAU^^ 
with a diagonal eigenvalue matrix A = diag(A/) and U G G1(A'^, C), then 

/(A) = [/diag(/(A,))[/-i. (1.4) 

If A cannot be diagonalized, a more general spectral definition can be derived from Eq. (1.3) using 
the Jordan decomposition [10, 1]. Non-Hermitian matrices typically have complex eigenvalues, 
and applying Eq. (1.4) to the sign function in Eq. (1.1) requires the evaluation of the sign of a 
complex number. The sign function needs to satisfy [sgn(z)]^ = 1 and, for real x, sgn(;c) = ±1 if 
X ^ 0. To satisfy these properties, it has become standard to define 

sgn(z) = 4t = sgn(Re(z)) > (1-5) 

where in the last equality the cut of the square root is chosen along the negative real axis. Using 
the definition (1.5) in the spectral definition (1.4) and reordering the eigenvalues according to the 
sign of their real part allows one to write the matrix sign function as 



sgn(A) = [/ ^ ,]U-'. (1.6) 
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The sign function satisfies sgn(A)^ = /, and a sliort calculation [7] shows that for this reason the 
overlap operator Dov(iU) as defined in Eq. (1.1) satisfies the Ginsparg-Wilson relation. Moreover, 
this definition agrees with the result obtained when deriving Eq. ( 1 . 1 ) from the domain- wall fermion 
formaUsm at / [11]. 

2. Arnold! method and function approximation for a non-Hermitian matrix 

A numerical implementation of the sign function using the spectral definition (1.4) is only 
possible for small matrices, as a full diagonalization becomes too expensive as the matrix grows. 
Alternatively, matrix-based iterative algorithms for the computation of the matrix sign function 
have been around for many years, see Ref. [12] and references therein. These are efficient for 
medium-sized problems, but are still unaffordable for the very large matrices occurring in typical 
lattice QCD simulations. Therefore, another iterative method is required which approximates the 
vector y = sgn(A)x, rather than the full sign matrix itself. Such iterative methods are already 
extensively used for Hermitian matrices [13, 14]. Most of these methods are derived from the 
Lanczos method, which uses short recuiTcnces to build an orthonormal basis in a Kiylov subspace. 

Krylov subspace methods have also been introduced for non-Hermitian matrices [15, 16]. The 
two most widely used methods to compute a basis for the Krylov subspace are the Arnoldi method 
and the two-sided Lanczos method. In contrast to the Hermitian case, the Arnoldi method requires 
long recurrences to construct an orthonormal basis for the Krylov subspace, while the two-sided 
Lanczos method uses two short recuiTcnce relations at the cost of losing orthogonality. Here we 
describe a Krylov subspace approximation based on the Arnoldi method to evaluate f{A)x for a 
generic function of a non-Hermitian matrix. 

We aim to construct an approximation to f{A)x using a polynomial of degree ^ — 1 with k <^N. 
For any k there exists a best polynomial approximation y = i {A)x of degree at most k—l, which 
is the orthogonal projection of f{A)x on the Krylov subspace Jtk{A,x) = span{x,Ax, . . . ,A'^^^x). 
An orthonormal basis Vi = (vi, . . . ,Vk) for the Krylov subspace Jik{A,x) is constructed using the 
Arnoldi recurrence 

AVk = VkHk + pkVk+iel , (2.1) 

where vi = x/ (5, (5 = \x\, Hj^ is an upper Hessenberg matrix, = and is the ^-th basis 

vector in C*^. Then VkVj^ is a projector on the Krylov subspace, and the projection y of f{A)x on 
J^k{A,x) can formally be written as 

y = VkV^f{A)x. (2.2) 

However, to compute the projection (2.2) one would already have to know the exact result f{A)x. 
Therefore, a method is needed to approximate the projected vector From Eq. (2.1) it follows that 

Hk = V^AVk , (2.3) 

which suggests the approximation [15] 

fm « V^f{A)Vk . (2.4) 
As ;c = pVke\, Eq. (2.4) can be substituted in Eq. (2.2), finally yielding the approximation 

y^l5Vkf{Hk)ei. (2.5) 
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In this approximation the computation of f{A) is replaced by that of f{Hk), where is of much 
smaller size than A. f{Hu)e\ should be evaluated by some suitable numerical method. 

The computation of the matrix sign function using Eq. (2.5) converges to the exact solution 
(see the m = curve in the left pane of Fig. 1). Unfortunately, in the case of the sign function, the 
convergence as a function of the size of the Kiylov subspace is very slow if some of the eigenvalues 
ai^e close to the function discontinuity along the imaginary axis. This problem can be resolved by 
deflation of these critical eigenvalues. 

For Hermitian matrices, it is well known that the computation of the sign function can be 
improved by deflating the eigenvalues smallest in absolute value [5]. Assume that m critical eigen- 
values Xi of A with orthonormal eigenvectors have been computed. Then 

m 

f{A)x = Y,f{Xi){u\x)ui + f{A)x^ , (2.6) 

!=1 

where x = x\\ +x± with = Y!i'Li{u]x)ui and x± = x — x^^. The first term on the right-hand side 
of Eq. (2.6) can be computed exactly, while the second term can be approximated using a Krylov 
subspace method for f{A)x±. Deflation will allow for a much smaller-sized Krylov subspace. 

For non-Hermitian matrices the eigenvectors are no longer orthogonal, and the simple decom- 
position into orthogonal subspaces, leading to Eq. (2.6), no longer holds. In the next two sections 
we will develop two alternative deflation schemes for the non-Hermitian case. 



3. Schur deflation 



We construct the subspace D.^ + J^k which is the sum of the subspace D.^ spanned by the 
right eigenvectors corresponding to m critical eigenvalues of A and the Krylov subspace ^(A,;c). 
Assume that m critical eigenvalues and right eigenvectors of A have been computed. From this, one 
can construct m Schur vectors Si, which form an orthonormal basis of Q.,,,, satisfying 

AS,n = S„jTf,j , (3.1) 

where Sm = {si, . . . ,Sm) and r„, is an m x m upper triangular matrix whose diagonal elements ai^e 
the eigenvalues corresponding to the Schur vectors. 

We propose a modified Arnoldi method to construct an orthogonal basis of the composite sub- 
space n„j + Jtk{A,x). That is, each Arnoldi vector is orthogonalized not only against the previous 
ones, but also against the Schur vectors si. In analogy to (2.1), this process can be summarized as 




siAvA 



A ( S„, Vk) = [ 5,„ Vt ) ( ""r," " 1 + PkVk+iel^k (3-2) 



Hk J 



with vi = x±/l5, where x± = {I — S,„Slj)x is the projection of x onto the orthogonal complement 
D.^ of n,„ and (5 = \x±\. The Hessenberg matrix 

^^(T,.SiAVk\ 
\ Hk 
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satisfies a relation similar to Eq. (2.3), namely H = Q^AQ, where the columns ofQ = {Sm Vu) form 
an orthonormal basis of n,„ + ^(A,x). In analogy to Sec. 2 we construct the approximation 

f{A)x^Qf{H)Q^x. (3.4) 

Because of the block structure (3.3) of H, the matrix f{H) can be written as 

where Y reflects the coupling between both subspaces and satisfies the Sylvester equation 

T^Y - YHu = f{T,„)X - Xf{Hk) (3.6) 
with X = SljAVk- Combining (3.4) and (3.5), we obtain 

f{A)x ^ S„J{T,„)Six+ {s^ Vk) [^f^fj^^ P^i ■ (3-7) 

f{T,„) and f{Hk)ei are computed with some suitable numerical method, and the mixed triangu- 
lar/Hessenberg Sylvester equation (3.6) for Y can be solved efficiently with the method of Ref. [1]. 



4. LR-deflation 

An alternative deflation in the same composite subspace Q.,n + ■y(fk{A,x) can be constructed 
using both the left and right eigenvectors corresponding to the critical eigenvalues. Assume that m 
critical eigenvalues of A and their con^esponding left and right eigenvectors have been computed, 

AR„, = R,„K,n , LlA = A„,4 , (4. 1) 

where A„, is the diagonal matrix of critical eigenvalues, and R,„ = (ri , . . . , r„,) and L,„ = (£i , . . . , ^„,) 
are the matrices with the corresponding right and left eigenvectors, respectively. The left and 
right eigenvectors corresponding to different eigenvalues are orthogonal. If the eigenvectors are 
normalized such that f-ri = 1, then L,^/?,„ = /,„, and /?,„L,^, is an oblique projector on the subspace 
Q.m spanned by the right eigenvectors. Let us now decompose xd&x = x\\+XQ, where x\\ = RmLj^x 
and Xq =x — x\\. Then 

f{A)x = R,„f{A,„)Lix + f{A)xe . (4.2) 

The first term on the right-hand side, which follows from Eq. (4.1), can be evaluated exactly, while 
the second term can be approximated by applying the Arnoldi method described in Sec. 2 to ;ce. An 
orthonormal basis Vjt is constructed in the Krylov subspace ^(A,;ce) using the Arnoldi recuiTcnce 
(2.1), with vi =Xq/P and (5 = \xq\. Successive operations of A on Xq will yield no contributions 
along the m critical eigendirections, hence ^(A,;ce) does not mix with D.,„. Applying the Arnoldi 
approximation (2.5) to Eq. (4.2) yields the final approximation 

f{A)x « R,nf{A,n)Llx + ^Vuf{Hu)ei . (4.3) 

Again, the first column of f{Hk) will be computed with some suitable numerical method. 
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Figure 1: Accuracy of approximation (4.3) for y = 5gn{H^{n))x with x = (1, 1, . . . , 1) at = 0.3, for a 4^ 
lattice (left) and a 6^* lattice (right). The relative error e = ||y — yll/Hyll is shown as a function of the Krylov 
subspace size k for various numbers of deflated eigenvalues m using the LR-deflation. 



5. Results 



We used the methods described above to compute the sign function occurring in the overlap 
Dirac operator ( 1 . 1 ) for a 4^^ and a 6"* lattice gauge configuration. Deflation of critical eigenvalues is 
essential because //w(m) has eigenvalues close to the imaginary axis. In practice, these eigenvectors 
need to be computed to high accuracy, as this will limit the overall accuracy of the function approx- 
imation. This was done with ARPACK [17]. The modified Arnoldi method was implemented in 
C++ using the optimized ATLAS BLAS library [18]. The convergence of the method is illustrated 
in Fig. 1 , where the accuracy of the approximation is shown as a function of the Krylov subspace 
size. The various curves correspond to different numbers of deflated eigenvalues, using the LR- 
deflation scheme. Without deflation (m = 0) the need for large Krylov subspaces would make the 
method unusable. Clearly, deflation highly improves the efficiency of the numerical method: as 
more eigenvalues are deflated, smaller Krylov subspaces are sufficient to achieve a given accu- 
racy. Furthermore, the deflation efficiency seems to grow with increasing lattice volume. Indeed, 
although the matrix size A'^ for the 6^ lattice is more than 5 times lai^ger than in the 4^ case, the 
Krylov subspace only has to be expanded by a factor of 1.2 to achieve a given accuracy of 10^^ 
(for m O.OOSA'^). It is also interesting to note that the modified Arnoldi approximation (4.3) for 
f{A)x is very close to the best approximation in the composite subspace, which is given by the 
orthogonal projection of f{A)x on Q.^ + J^k{A-,x), as was checked numerically. 

The results for the Schur deflation are not shown here, but are very similar to those for the LR- 
deflation. The Schur deflation is slightly less accurate, and requires more CPU time per evaluation, 
mainly because of the additional orthogonalization of the Arnoldi vectors with respect to the Schur 
vectors. However, the time taken by its initialization phase is halved, as it only requires the compu- 
tation of the right eigenvectors, and the best choice of deflation scheme will depend on the number 
of vectors x for which sgn(//w)-^ needs to be computed. If one needs to apply both sgn(//w) and its 
adjoint, then, obviously, the LR-deflation will be the better choice. A more detailed discussion of 
both deflation schemes can be found in Ref. [1] 
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6. Conclusion 

In this talk we presented an algorithm to approximate the action of a function of a non- 
Hermitian matrix on an arbitrary vector, when some of the eigenvalues of the matrix lie in a region 
of the complex plane close to a discontinuity of the function. The method approximates the solution 
vector in a composite subspace consisting of a Krylov subspace augmented by the eigenvectors cor- 
responding to a small number of critical eigenvalues. Two deflation variants were presented based 
on different subspace decompositions: the Schur deflation uses two coupled orthogonal subspaces, 
while the LR-deflation uses two decoupled but non-orthogonal subspaces. Deflation explicitly 
takes into account the contribution of the critical eigenvalues. This allows for smaller-sized Krylov 
subspaces, which is crucial for the efficiency of the method. The method was appUed to the overlap 
Dirac operator of lattice QCD at nonzero chemical potential, where the importance of deflation was 
clearly demonstrated. 
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